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Abstract 

The dynamics of inelastic hard spheres is described in terms of the binary 
collision expansion, yielding the corresponding pseudo-Liouville equation and 
BBGKY hierarchy for the reduced distribution functions. Based on cluster 
expansion techniques we derive the Boltzmann and ring kinetic equations for 
inelastic hard spheres. In the simple ring approximation, we calculate the 
structure factor of vorticity fluctuations in a freely evolving, dilute granular 
gas. The kinetic theory result agrees with the result, derived previously from 
fluctuating hydrodynamics. In the limit of incompressible flow, this structure 
factor alone determines the spatial velocity correlations, which are of dynamic 
origin and include long range r _ci -behavior. The analytic results are compared 
with MD simulations. 

I. INTRODUCTION 

In recent years the interest in static and rheological properties of assemblies of mesoscopic 
or macroscopic bodies or granules has been rapidly increasing. Both as a solid and as 
a fluid, granular systems show unusual behavior []!]. Flows of granular material can be 
subdivided into quasi static (contact) flows and rapid (collision driven) granular flows. In 
this characterization of Ref. fl|], the former is referred to as the granular liquid regime and 
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the latter as the granular gas regime. Only in rapid granular flows, the dynamics can be 
described by sequences of binary collisions and, as a consequence, the methods of kinetic 
theory are most suitable p|-|TTp. Such flows do obey the standard conservation laws of 
mass and momentum, and can therefore be considered as fluids. However, energy is not 
conserved. The dynamics is essentially dissipative, which gives rise to several interesting 
new phenomena, such as clustering and inelastic collapse U. 

In rapid granular flows, collisions of granules are accompanied by conversion of kinetic 
energy into rotational energy and into energy of other internal degrees of freedom, and lead 
to effective 'cooling' phenomena in unforced flows. Such flows have been extensively modeled 
through smooth and rough hard spheres with inelastic collisions pHTl[ • 



Up to now the theoretical description of rapid granular flows has been based on 
Boltzmann-Enskog kinetic equations for hard sphere-type interactions. Inherent to such 
descriptions is the molecular chaos assumption of uncorrelated binary collisions. To study 
the extent to which collective effects may be of importance at a fundamental level of de- 
scription, one needs to correct for the breakdown of the molecular chaos assumption and to 
account for the effects of dynamic correlations. 

In the last 35 years many-body theories have been developed to account for these dynamic 
correlations in systems of microscopic particles obeying the standard conservation laws. 
The fundamental concept to describe these dynamic correlations are 'ring collisions', i.e. 
sequences of correlated binary collisions, which lead to the so called ring kinetic theory. 
This ring kinetic theory for systems of perfectly smooth elastic hard spheres has been at 
the basis of all major developments in nonequilibrium statistical mechanics over the last 
three decennia: it explains the logarithmic density dependence of the transport coefficients, 



and the breakdown of the virial expansion for transport coefficients |L2|-|l4j]; it explains the 
algebraic long time tails of the velocity autocorrelation function and similar current-current 
correlation functions [p)|-|T7|j , the non-analytic dispersion relations for sound propagation 



and for relaxation of hydrodynamic excitations [|18|| , as well as the breakdown of the Navier- 
Stokes equations in two-dimensional fluids at very long times and the non-existence of linear 



transport coefficients in 2-D ||19|| ; moreover, it explains the existence of long range spatial 
correlations in nonequilibrium stationary states p0|j2l[] , driven by reservoirs which impose 
shear rates or temperature gradients, or in driven diffusive systems. Such systems violate 



the conditions of detailed balance and the stationary states are non-Gibbsian states |22 |. 

The goal of this paper is to explain the long range correlations (in fact, intermediate range, 
as the algebraic tails are exponentially cut off at the largest scales), observed in molecular 
dynamics simulations of unforced flows of two-dimensional granular gases |1TT1J23| , |2^| . This 
will be done by extending the mean field-type Boltzmann-Enskog equation for rapid granular 
flows by including ring collisions. To do so we choose the model of perfectly smooth, but 
inelastic hard disks or spheres (d = 2,3) of diameter a. The model incorporates the most 
fundamental feature of the dissipative dynamics of granular flow, namely the conversion of 
kinetic energy into internal energy. Moreover, it is a many-body system with well-defined 
and relatively simple dynamics, which has already widely been used in molecular dynamics 
simulations [0-0,0,0. 

The interactions between smooth inelastic hard spheres (IHS) only affect the translational 
degrees of freedom and are modeled by instantaneous collisions as in the case of elastic hard 
spheres. During a collision momentum will be transferred along the line joining the centers 
of mass of the two colliding particles, indicated by the vector er pointing from the center of 
particle 2 to that of particle 1, as illustrated in Fig. 1. 

Upon collision, inelastic hard spheres of equal mass m change their velocities according 
to the collision rules 

v* = vi - |(1 + a)(vi2 • <t)<t 

v 2 * = v 2 + i(l + a)(vi 2 -<x)<x, (1) 

where er = cr/a denotes a unit vector, v 12 = Vj — v 2 , and a is the coefficient of normal 
restitution. Whereas the total momentum of the two particles is conserved in a collision, 
the total energy loss is |me(vi2 • cr) 2 , where the coefficient of inelasticity is e = 1 — a 2 . 
The restituting (precollision) velocities (v**,v|*) giving rise to (v 1; v 2 ), are found by 



inverting collision rule ([I]) and given by (see Fig. 1) 

vr = vi-|(l + a- 1 )(v 12 - < r)<T 

v 2 ** = v 2 + l(l + a' 1 )(v 12 -&)&. (2) 

Note that this inversion is not possible if a = 0, which is a trivial limit. 

The kinetic theory for inelastic hard spheres will be developed in close analogy with that 
for elastic ones. The rather singular streaming operators St(T), which generate the T-space 
trajectories of the iV-hard sphere system, are reformulated (see section [□]) in terms of T- 
operators (binary collision operators), following the original derivation in Ref. |25). This 
allows us to introduce a pseudo-Liouville equation for the iV-particle T-space distribution 
function, and obtain the Bogoliubov, Born, Green, Kirkwood and Yvon (BBGKY) hierar- 
chy for the reduced distribution functions. The Liouville equation or the BBGKY hierarchy 
forms the standard starting point for deriving kinetic equations, such as the Boltzmann equa- 
tion and the ring kinetic equation. This program is carried out in section [[II]. Subsequently, 
we use the ring kinetic equation for inelastic hard spheres to calculate the nonequilibrium 
pair distribution function in a freely evolving, dilute granular gas. 

A particular solution of the Boltzmann equation for this system is the so called homo- 
geneous cooling state (HCS), where the velocity field vanishes everywhere, and the density 
n and temperature T(t) are spatially homogeneous, while the temperature decays in time. 
Some basic features of this HCS-solution are summarized in section |V] They are necessary 
for an understanding of the remaining part of this paper, including its instability against 
spatial fluctuations in flow velocity and density. In section |V| we will calculate the build- 
up of vorticity correlations using the ring kinetic equation. The resulting expressions for 
the spatial correlations (u a (r + r', t)up(r', t)) in the flow field, including long range r~ d - 
behavior, are identical to the results from fluctuating hydrodynamics for incompressible 
velocity fluctuations |23| in the low density limit. The incompressibility assumption gives a 
valid description of the velocity correlations up to a length scale that diverges as 1/e, beyond 
which the r _d -tail is cut off exponentially. Section [V| concludes with a comparison of the 
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theoretical predictions with the results of molecular dynamics simulations of unforced flows 
in a two-dimensional gas of inelastic hard disks (see |?|f||24]|). 



II. THE BINARY COLLISION EXPANSION 

Having defined the dynamics of the system, we now turn to the statistical mechanics. 
The ensemble average of a dynamical quantity A(T) obeys the equality 

J dTp(T, 0)A(V(t)) = J dTA(T)p(T, t), (3) 

where T = {xi,x 2 , ■ ■ -x^} with x« = {r-j, v.;} is a point in the iV-particle phase space. On 
the left hand side the time dependence is assigned to the dynamical variable A(T(t)) = 
S t (X)A(T) with S t {T) the time evolution operator, and on the right hand side to the N- 
particle distribution function p(T,t), which can be done by considering the terms in (|3|) as 
an inner product. The time evolution of the iV-particle distribution function is then given 

by 

p(T,t) = Slp(T,0), (4) 

where Sj(T) is the adjoint of S t (T). For Hamiltonian dynamics Sj(T) = S- t (T). 

For particles with hard core interactions the dynamics is undefined for physically in- 
accessible configurations, where the particles are overlapping. Such configurations have 
a vanishing weight in Eq. (|3p. Since St(T) only appears in the combination p(T,0)St(T) 
which vanishes for overlapping initial configurations, it suffices to consider W^(T)S t (T), 
with ^^(r) = n«<j Wfaj) where W(r) is the overlap function: 

!0 if r < a (overlapping) 
(5) 
1 if r > a (non-overlapping). 

However, the methods of many-body theory require formal perturbation expansions and 
subsequent resummations. To do so, the time evolution operator S t (T) needs to be defined 
for all configurations, including the unphysical overlapping configurations. A convenient 



representation, defined in all points in phase space, has been developed for elastic hard 



spheres in Ref. p5], and is based on the binary collision expansion of S t (T) in terms of 
binary collision operators T{ij). 

The binary collision operator may be defined in terms of two-body dynamics through 



the time displacement operator #(12), given in Ref. by 

#(12) = fi?(12) + f dr# (12)T(12)#°_ r (12). (6) 

J 

Here the free streaming operator S?(12) = exp[tL°(12)], with L°(12) = L\ + L\ and L° = 
Vj ■ d/dri. Its action on an arbitrary dynamical function is 

S?(12)A(r x , p 1; r 2 , p 2 ) = A(r a + v^, Pl , r 2 + v 2 t, p 2 ). (7) 



Following the argument of Ref. [25| for the case of elastic hard spheres step by step, the 



binary collision operator T(12) for inelastic hard spheres is constructed as 

T(12) = o d - x f d<r|v 12 ■ &\8{r 12 - - 1), (8) 



' V12-CT<0 

where 6* is an operator that replaces all (precollision) velocities Vj (i = 1,2) appearing to 
its right by postcollision velocities v*, defined for the inelastic case through collision rule 
(HI), and d is the dimensionality of the system. 

The binary collision operator is defined for overlapping and non-overlapping configura- 
tions of two hard spheres. It extends the definition of #(12) to all points in phase space. 
In the ensemble average considered in Eq. (^), the overlap function Wn(T) contains a fac- 
tor W(r 12) which vanishes whenever r 12 < a. Moreover, the pseudo-dynamics introduced 



through Eq. @ is noninvertible (see [p5] ). Consequently, the generator #(12) for two- 
particle dynamics is only defined for positive times. For later reference we quote the prop- 
erty T(12)v4(12) = 0, where A(12) is a constant or a function of the argument (v x + v 2 ), 
because of conservation of particle number and linear momentum in binary collisions. As 
kinetic energy is not conserved in inelastic collisions, this property does not apply if A(12) 
is a function of the argument (vf + u|). 
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Combinations of T-operators and free streaming operators S®, preceded by appropriate 
combinations of overlap functions, can be used to construct the time displacement operators 
St for dynamical variables. However to discuss the Liouville equation and describe the 
time evolution of the reduced distribution functions we need to consider the adjoint time 
displacement operators. 

We start with the two-particle streaming operator St(12). In order to obtain the adjoint 
S} (12), we consider the integral equality 

J dxxdx 2 B(12)W(12)S t (12)A(12) = J dx 1 dx 2 A(12)S 1 t (12)W(12)B(12), (9) 

where dxidx 2 — dr 1 dv 1 dr 2 dv 2 and A(12) and -8(12) are arbitrary functions of the phases 
x\ and x 2 . Since VK(12) is appearing to the left of St(12), this integral is well-defined. 
Substituting St(12) from Eq. (|j) and using Liouville's theorem S%(12) = S®(12) for free 
particle motion, the left hand side of Eq. can be written as 

J d Xl dx 2 {A(12)S°_ t (12)W(12)B(12) + 

f dr[S°_ T (12)W(12)B(12)]T(12)Sl T (12)A(12)}. (10) 

Now notice that the binary collision operator T(12) contains two terms: a real collision 
term with b* and a virtual collision term. The real collision term appearing in fllOl ) can 
be transformed using v* 2 ■ & = -av 12 ■ <r, yielding a Jacobian dvxdv 2 = dv^dv^/ct, and 
then relabeling v* — > Vj and Vj — > v**. In the virtual collision term we relabel & — > —a 
and use the free particle Liouville theorem to finally write (|ToD in terms of the adjoint 
time displacement operator, defined through the right hand side of Eq. @. The resulting 
expression for the adjoint time displacement operator becomes 

4(12) = 5^(12) + ^^5^(12)^(12)^^,(12), (11) 

where the binary collision operator T is defined as 

T(12) = a d - 1 [ d<r(v 12 • a) (^S(r 12 - - 5{r 12 + a)) . (12) 
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Here &** acts on the velocities v, (i = 1, 2) to its right and replaces them by restituting ones, 
v**, as defined in collision rule (0). 

A property of T-operators, equivalent to T(12)v4(12) = 0, is 

J dxidx 2 A(12)T(12) J B(12) = 0, (13) 

if A(12) is a constant or a function of the argument (v x + v 2 ), whereas -8(12) is an arbitrary 
function. 

The time displacement operators St(12) can be put in a more convenient form by using 
the property, T(12)S' t (12)T(12) = 0, valid for any t > 0. It also holds with T replaced by T. 
This relation expresses the fact that two hard spheres cannot collide more than once with 
only free propagation in between. Using this property, the time displacement operators can 
be written as 

W(12)S t (12) = W(12) exp[tT°(12) + tT (12)] 
Sf(12)W(12) = exp[-tL°(l, 2) + tT(12)]W(12). (14) 

We now return to the full iV-particle system. As shown in Ref . , the dynamics of the 
iV-particle system can be represented in the compact form of pseudo-streaming operators 
with the help of the above property T(12)S' i (12)T(12) = 0, yielding 

W N (T)S t (T) = W^(r)exp[tL°(r) +tY J T{ij)} 

i<j 

Sj(T)W N (T) = e W [-tL (T)+tJ2TmW N {r), (15) 

i<j 

with L°(T) = J2iL® the free particle streaming operator. These time evolution operators 
are defined everywhere in phase space, and the overlap function gives a vanishing weight to 
unphysical configurations, provided that Wn(T) always appears to the left of T-operators, 
or to the right of T-operators. 

A minor generalization is to include a conservative external force field in the dynamics. 
In that case the single-particle free streaming operator should be defined as 

^ = ^ — + 8,. A (16) 
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where is the external force per unit mass, acting on the z-th particle. 

Similar results for the binary collision operators T and T for inelastic hard spheres have 
been derived independently by Brey et al. 0. 

III. BBGKY HIERARCHY 

The time evolution of the iV-particle distribution function p(T,t) = p(x±, x 2 , ■ ■ ■ x^, t), 
with Xi = {rj, Vj} is given by the Frobenius-Perron equation. For conservative Hamiltonian 
systems this equation is the Liouville equation which is an expression of the incompressibility 
of the flow. In the case of dissipative systems which by definition are time irreversible, the 
phase space volumes are contracted along the flow. According to Eqs. (§) and ([15|) the 
time evolution of the distribution function for inelastic hard spheres is given by the pseudo- 
Liouville equation 

[d t + L°(T)}p(r,t) = j:m)p(r,t). (17) 

i<j 

An equivalent representation of the time evolution of the system can be given in terms of 
reduced s-particle distribution functions (s = 1,2,...), defined as 

N\ r 

/i2... s (*) = P s) (x 1 ,x 2 , ...x s ,t)= — - J dx s ... dx N p(x u x 2 ,... x N , t), (18) 

where p(t) is normalized to unity. Integration of Eq. (0) using fllTf ) yields the BBGKY 
hierarchy for the reduced distribution functions. We only quote the first two hierarchy 
equations: 

{d t + Ll)h = J dx 2 T(12)/ 12 
[d t + Li + L° 2 - T(12)]/ 12 = J dx 3 [T(13) + T(23)]/ 123 . (19) 

This set of equations is an open hierarchy, which expresses the time evolution of the s-particle 
distribution function in terms of the (s + l)-th function. 

In the literature on kinetic theory of inelastic hard spheres the first equation of the 
BBGKY hierarchy has frequently been derived intuitively and used as a starting point to 
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obtain the Boltzmann-Enskog equation for the single-particle distribution function [Q,[||5]||. 



Using the explicit expression p2|) for T(12) it can be written in full detail as 



er x 



(d t + L?)/(r X) vi, t) = o d ~ l f dv 2 / d<x(v 12 
{^/ (2) (ri, < , ri - <x, vT, f) - / (2) (n, Vi, n + <r, v 2 , f)} . (20) 

The corresponding equation for the rate of change of an average is obtained by multi- 
plying the first hierarchy equation in (|19|) with J dv 1 , 0(v 1 ), and using the adjoint of T(12) 
to find 

d f d f 

— J dv 1 ^(v 1 )/(r 1 , vi, t) + — • J dv 1 v 1 ^(v 1 )/(r 1 , v l5 1) 
= y dvi I dx 2 / 12 T(12)V(v 1 ) 
= / dvi / dv 2 / d<r(v 12 • <x) x 

/( 2 )(r 1; vi, ri + a, v 2 , t)[^K) - ^( Vl )], (21) 

where we have set the external force equal to zero. This equation is the starting point for 
deriving macroscopic conservation laws, hydrodynamic equations and the rate of change of 
the temperature. 

We return again to the kinetic equations. In order to derive a closed equation for the 
single-particle distribution function / or for the pair function f\ 2 some kind of closure relation 
is required, to express fu in terms of /, such as Boltzmann's molecular chaos assumption, 
Bogoliubov's functional assumption [p^] or cluster expansion methods Here we shall 



illustrate how the methods derived in the kinetic theory for elastic hard spheres can be 
transferred directly to the case of inelastic hard spheres. This will be done by deriving the 
Boltzmann equation and the ring kinetic equation for the model under consideration. 

The Boltzmann equation is obtained from the first hierarchy equation by keeping only 
terms to dominant order in the density. This implies the fundamental assumption of molec- 
ular chaos for dilute gases, expressing the absence of dynamic correlations between the 
velocities of two colliding particles just before collision, i.e. at | r 12 1 = a + 0, 
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/i2 = /(r 1 ,vi,t)/(r 2 ,v 2 ,t). (22) 

Furthermore, in the low density limit the spatial separation between the colliding particles 
can be neglected, and the binary collision operator T(12), entering in the BBGKY hierarchy, 
reduces to 

T(12) = 5(r 12 )T (12) = 8(r 12 ) a^ 1 f d<x(v 12 ■ &) (\b* a * - l) . (23) 
Then the nonlinear Boltzmann equation for inelastic hard spheres becomes 
{d t + L° l )f(r 1 ,v 1 ,t) = a d - 1 ( dv 2 / d<r(v 12 • &) x 

{^/(ri, vr , t)/(ri, vT, t) - /(n, vx, t)f(r u v 2 , t)} = /(/, /). (24) 

There are two significant differences with the Boltzmann equation for the elastic case: (i) 
the occurrence of 1/a 2 in the gain term on the right hand side of (|24]); one factor 1/a 



comes from the Jacobian dv^dvg* = (l/a)dvidv 2 and the other one from the reflection law 
• cr = — (l/a)v 12 • a (see Fig. 1). (ii) In the inelastic case, the restituting precollision 
velocities, which yield (v 1; v 2 ) as postcollision velocities, are different from the postcollision 
velocities (v*,v 2 ), which result from the direct precollision velocities. In the elastic case 
(a = 1) the relation v* = v** holds. 

The Boltzmann- Enskog equation for inelastic hard spheres (see Refs. [@,|],|5]|1) is obtained 
by replacing f 12 in the first hierarchy equation of ( |i~9l) by x(r l5 r 2 |n)/!/ 2 , where x(r 1; r 2 |n) 
is the pair correlation function of elastic hard spheres in a spatially nonuniform equilibrium 



state (see Ref. [0). This version of the molecular chaos assumption still neglects the 
velocity correlations, built up by sequences of correlated binary collisions, but does account 
for static short range correlations, caused by excluded volume effects. The consequence for 
the transport coefficients have been worked out in Refs. 

As the density increases the contributions of correlated collision sequences to the col- 



lision term on the right hand side of ( |20| ) become more and more important. The 



most simple sequence of correlated collisions are the so called ring collisions; for exam- 
ple (12)(13)(14). . ,(23')(24'). . .(12), ending with a recollision of the pair (12), which was 
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involved in the first collision. In the intermediate time particle 1 collides, say, s times with 
s different particles (3,4,. . .), and particle 2 collides s' times with another set of s' different 
particles (3', 4',. . .). When particles 1 and 2 are about to recollide, they are dynamically 
correlated through their collision history, and the molecular chaos assumption (^) is no 
longer valid, i.e. g 12 = fi 2 - /1/2 0. 

A simple way to take these correlations into account at moderate densities has been given 
in Refs. p7"y2"9fl . The method is based on a cluster expansion of the s-particle distribution 
functions, defined recursively as 

/12 = /1/2 + #12 

/l23 = /1/2/3 + f 19-23 + /2.913 + /3.912 + 0123, (25) 

etc. Here gu accounts for pair correlations, #123 for triplet correlations, etc. The molecular 
chaos assumption implies g\ 2 = 0, which is equivalent to fl22|) . The basic assumption to 
obtain the ring kinetic equation is that the pair correlations are dominant and higher order 
ones can be neglected, i.e. gi 23 = (71234 = • • • = in cluster expansion fl2"5]). 

Substitution of ( P5] ) into (|T^) and elimination of dfi/dt (i = 1,2) from the second hier- 
archy equation using the first one, yields the ring kinetic theory of inelastic hard spheres: 



(d t + L^h = J dx 2 T(12)(f 1 f 2 + g 12 ) (26) 
[d t + L\ + L G 2 - T(12) - (1 + V n ) J da; 3 T(13)(l + V Vi )h]g l2 = T(12) [fj 2 + g 12 ) . (27) 

Here V{j is a permutation operator that interchanges the particle labels i and j. The second 
equation is the so called repeated ring equation for the pair correlation function. If the 
operator T(12) on the left hand side of the second equation is deleted, one obtains the 
simple ring approximation. Formally solving this equation for g 12 yields an expression in 
terms of the single-particle distribution functions fi (i — 1, 2, 3), and subsequent substitution 
into the first hierarchy equation above yields the generalized Boltzmann equation in ring 
approximation. For a more detailed discussion of the collision sequences taken into account 
by Eq. (|2~7D we refer to the original literature [p7lp9l, [3~0|. 
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The kinetic equations (|26| ) and (|27|) constitute the new extensions of this article: the ring 
kinetic equations for inelastic hard spheres. All standard results for elastic hard spheres are 
recovered by setting the restitution coefficient a = 1. In this paper we will determine from 
(p7|) the behavior of the nonequilibrium pair distribution function g(xi,X2,t) on hydrody- 
namic time scales, for a special solution f(x,t) of the Boltzmann equation. The spatial 
correlations (u a (r + r', t)u/3(r', t)) of the flow field u(r,t), calculation of which is our direct 
goal in this paper, can be directly obtained from g(x\,X2,t) by integration. In the next 
section will summarize some known results about the homogeneous cooling state which are 
necessary for an understanding of what follows later on. 

IV. HOMOGENEOUS COOLING STATE 

A gas of elastic hard spheres will relax to local equilibrium on a (kinetic) time scale of a 
few mean free times, t ~ Io/vq, where Iq is the mean free path and Vq the thermal velocity, 
i?m>Q = T. Here the Boltzmann constant ks is set equal to unity. Finally, the system will 
reach a spatially homogeneous equilibrium state on a (hydrodynamic) time scale ~ L/vo, 
where L is the linear dimension of the system. 

However, in the case of inelastic hard spheres kinetic energy is lost in collisions, and if 
the system is not driven, the kinetic energy associated with the thermal motion decreases, 
and interesting instabilities occur, such as clustering |7j and inelastic collapse 0. 

Here we are interested in a special solution of the Boltzmann equation ([24|) for inelastic 
hard spheres, the so called homogeneous cooling state (HCS). In this state the distribu- 
tion function /(r,v,t) = f(y,t), as well as the hydrodynamic fields are spatially homoge- 
neous. These functions are defined as density n(r,t) = J dv/(r, v, t), flow velocity u(r, t) = 
(l/n(r,t)) J dvv/(r, v, t), and granular temperature T(r,t) = (m/dn(r,t)) J dvK 2 /(r, v, t), 
where V(r, t) = v — u(r,t) is the peculiar velocity. Furthermore, the flow velocity can 
be taken to vanish, u(r, t) = 0, and n(r, t) — n is constant in space and time. However, 
T(r, t) = T(t) depends on time. Based on the fundamental concepts of the Chapman- Enskog 



13 



theory, we expect that the single-particle distribution function / after an initial transient of 
the order of a few mean free times, will only depend on time through its first few moments, 
which is here only the temperature T(t). For dimensional reasons /(v, t) then takes the 
scaling form 

f ^=wA^)- (28) 

with the thermal velocity Vo(t) depending on time. One can derive an integral equation 
for the unknown scaling form /(c), with c = v/vo(t), by inserting fl28f ) into the Boltzmann 



equation (|24[) . The result is 



-^(" + e >-^)^)=-»*-W,/j, (29) 

with 

/(/, /) = /dc 2 / d<r(c 12 ■ Or) (-J™ - l) /( Cl )/(C 2 ). (30) 

J Jc 12 -a>o Vcr J 

To determine the rate of change of the temperature, we use Eq. fl2T|) with ^(v) = v 2 and 
fi2 = f( v i,t)f(v2,t), and calculate v\ 2 — v\. With the help of Eqs. ([I]) and (P^), we obtain 

% = -^mno^vfo = -2^ 7 T, (31) 

where u is the time dependent collision frequency given in Boltzmann theory by uo = 



fidW d ~V T / 7rm with tt d = 2n d / 2 /T(d/2) the surface area of a d- dimensional unit sphere, 
and 7 is a time independent cooling rate, defined by 



= iS) ^/ dCl / dC2 i,, S >o d * (Cl2 '* )3/(C ' )/(C2) ' (32) 
which depends on the unknown scaling form /(c). 

At this stage, it is convenient to change to a new time variable r, defined through 
dr = u(T(t))dt. Integration yields 

T = -ln(l + 7t/*o)- (33) 
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The time r actually presents the average number of collisions suffered per particle. In the 
elastic limit (e — > 0), it becomes the real time, r = t/t , measured in units of the mean free 
time to = 1/^(^0 ) at the initial temperature T . 

Eq. ( |3T| ) is readily integrated to obtain for the temperature 

r| " = (uW =T,eip( " 2ir) ' (34) 

The above equation represents the well known algebraic decay law for the granular temper- 
ature in the homogeneous cooling stateQ (see Refs. |7]|| and references therein). 

Combining ( f29|) and (|3lD yields an integral equation for the scaling form /(c), i.e. 

_^ +Cl ._y /(ClH ?(/,/,, ( 35) 

where 7 is given through (|32l ) in terms of / itself. By a moment expansion of this equation, 



it has been shown |3l] that the fourth cumulant of / is small for any value of the inelasticity 
e both in d = 2 and 3. This has been confirmed by a direct simulation Monte Carlo calcu- 
lation of the fourth and sixth cumulant in the homogeneous cooling state of inelastic hard 



spheres (d = 3) |j32[ . To a good approximation, therefore, the solution of fl35|) approaches 
a Maxwellian, i.e. f(c) ~ 0(c) = n~ d / 2 exp(— c 2 ). Similarly, the dimensionless cooling rate 
(j32|) is well approximated by replacing /(c) by 0(c) in 



^—^h^^^-if- < 36) 



V. RING KINETIC THEORY 

In this section we will derive the structure factor of transverse velocity or vorticity 
fluctuations, S±(k,t), for a freely evolving, dilute gas of inelastic hard spheres. In general 



1 Note that, because of the spatial homogeneity of /, this result is still valid in Enskog theory, for 
general densities, once the collision frequency has been increased by the factor x( n )> which is the 
local equilibrium pair correlation function at contact. 
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velocity correlations can be described by an isotropic tensor S a p(\i,t), related to the scalar 
isotropic functions S±(k,t) and Su(k,t) via the decomposition 



S a/3 (k,t) = k a kpS\\(k,t) + (5 a p - k a kp)Sx(k,t), (37) 

where hats denote unit vectors. The structure factors S±(k, t) and S\\(k, t) of transverse and 
longitudinal velocity fluctuations are also related to the energy spectrum function E(k) in 
the theory of homogeneous turbulence |33] . 



The inverse Fourier transform of S a p(h, t) are the spatial correlations G a p(r, t). In terms 
of the microscopic velocity field u(r, t), it may be defined as 

G a p{r, t) = -J dr'( MQ (r + r', t)up(r', t)). (38) 



A similar decomposition holds for the spatial velocity correlations: 

G a p(r,t) = r a rpG\\{r,t) + (S a p - r a f-p)G ± (r,t). (39) 

Whereas the tensors G a p(r,t) and S a p(k,t) include self-correlations of particles, the corre- 
lation functions gi2 and S12 (as defined below) occurring in the second hierarchy equation 
do not. Therefore, it is convenient to substract the self-correlation part and introduce the 
tensors S^p(k,t) = S a p(k,t) — T(t)5 a p/nm and C^g(r, t) = G a p{v,t) — T(t)5 a p5(r)/nm, 
which is regular at the origin and related to by 

Gtp( r , t ) = yJ dr '^ / dvi / d ^ViaV 2 p9{r + r', v 1; r', v 2 , t). (40) 

Incompressibility of velocity fluctuations then implies S^(k,t) = 0, in which case G\_{r, t) is 
related to G^(r,t) by (see f33[) 

Gi(r,t) = q(r,t) + ^ I |-q(r,t). (41) 

As discussed in section |ITj], in the low density limit, T(12) = 5(ri 2 )T (12), the first 
hierarchy equation reduces to the Boltzmann equation, and the second hierarchy equation 
to its simple ring approximation, 
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di + Vl 'fr~ 1 +Ql + Y2 'fr~ 2 +n2 ) 9l2 = 5 ( r i2)^o(12)/i/ 2 , (42) 
where Q is minus the linearized Boltzmann collision operator, i.e. (i = 1, 2) 



-a*" 1 / dv 3 / d<x(v i3 • or) (- 2 b*; - l) X 

{/(v J ,t)^(v J ) + /(v3,^(v 3 )}. (43) 



The term on the right hand side of Eq. ([42]), which vanishes in the case of detailed balance, 
provides a source of correlations for inelastic hard spheres in the homogeneous cooling state. 
For a homogeneous distribution /(v, t), g(ri, vi, r 2 , v 2 , t) depends on r i2 = ri — r 2 only, and 
the ring equation for the Fourier transform s(k, Vi,v 2 ,t) = / dri 2 e _lk ' ri2 g(ri, vi, r 2 , v 2 , t) is 
given by 

(d t + ik • v 12 + ^ + fi 2 )s 12 = T„(12)/(vi, t)/(v 2 , t). (44) 



Inspired by the scaling ansatz ( pB"D for the homogeneous cooling solution of the Boltzmann 
equation, we write the density and temperature dependence of s i2 explicitly as 

s(k, vi, v 2 , t) = s(k, ci, c 2 , r). (45) 

Again it is convenient to transform to the kinetic time r, so that Eq. ([44]) reduces to 

[d T + Aia(k)]3(k, c x , c 2 , r) = / a^ 1 r (12)/(c 1 )/(c 2 ), (46) 

with the mean free path Z — ^o/ 1 ^- Here we have introduced the notation (j = 1, 2) 



Aj-(k) = il k ■ cj + 7 ( d + Cj ■ + l^na^Q, 



dc o. 



A 12 (k)=A 1 (k)+A 2 (-k), (47) 



defined in terms of the dimensionless quantities 



T (12) = 1 ,,^ 0(12). (48) 
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The formal solution for the case that all pair correlations vanish at the initial time t = 0, 
is given by 

S ~(k, Cl ,c 2 ,r) = 1 " 6X P [ ~;f 12(k)] Z ^- 1 r (12)/(c 1 )/(c 2 ) 

Aia(k) 

a,/x zaV«-j -+- z^k; 

(^ A L (c 1 ,k)^(c 2 ,k)|/ 0( T d - 1 f (12)/(c 1 )/(c 2 )). (49) 

In the second line we have made a decomposition in eigenfunctions of the operator A and 
used the bracket notation, with inner products representing integrals over c. The left and 
right eigenfunctions satisfy the eigenvalue relations, 

A,(k)|^(c„k)) = -z A (k)|^( Cj ,k)) 

(^(c„k)|A,(k) = -z x (k)^ L ( Cj ,k)\- (50) 
The structure factor of transverse velocity fluctuations, S±(k,t), is then given by 

S±(k,t) = ~2" / dv l / dv 2^1±a^2± a s(k, Vi, V 2 ,t) 

= Vo(t){ci± a C2± a \s(k,Ci,c 2 ,r)). (51) 

We state here without further derivation that the transverse velocity or shear modes are the 
only modes contributing to S±(k, t) (see also pT| , p3| ). It can be shown that the eigenfunctions 



corresponding to the shear mode are ip± a = c ■ kj_ a = c± a , where the subscript a denotes 
one of the d — 1 degenerate shear modes. The dispersion relation for the relaxation rate 
z±(k) of the shear mode with wavenumber k is z±(k) = 70(1 — k 2 ^ 2 ). The correlation 



length £ = Jrj/pu jo with rj ~ yT(t) the time dependent shear viscosity, p = mn the 



mass density, lu ~ yT(t) the collision frequency and 70 = e/2d. For small inelasticity 
£ diverges as l/y/e. So in calculating S±(k,t), we only have to calculate the quantity 
(ci± a c 2 _ La |/ o" a! ~ 1 To(12)/(ci)/(c 2 )) = -Q d l cr d - 1 'y/y/2ir = -7/n « -70/n, as can easily be 
shown using Eq. (|2~T|). We then obtain for S^(k,t) the expression 

S|(M) = (iw) !^i(iz|aw . (52) 

V nm J 1 — 



This low density result agrees with a previous result derived from fluctuating hydrodynamics 
|23| , which, moreover, extends the validity of the above expression to higher densities. 

The subsequent analysis of G±(r, t) and G»(r,t), under the simplifying assumption that 
the fluctuating flow fields are incompressible, i.e. St(k,t) = 0, has been given in P3"| . The 
full analysis for the compressible case can also be given, and yields only small modifications 
to the incompressible case, except at the largest scales, where the algebraic tails are cut off 
exponentially |3l| . 



It is found that the spatial correlation functions G±(r,t) and G\\(r,t) show structure 
on spatial scales large compared to the mean free path Iq. In particular one obtains for 
asymptotically large r the algebraic tails, G»(r, t) ~ — (d — l)G±(r, t) ~ A/r d with explicit 



expressions for A |3T|. Furthermore, G»(r, t) is positive everywhere, while G±(r,t) has a 
negative minimum, reflecting the presence of vortices in the system. The structure function 
S±(k, t) for the flow field in the inelastic hard sphere system is the analog of the energy spec- 
trum function E(k) in the theory of 2-D or 3-D homogeneous turbulence in incompressible 
fluids. The qualitative forms of G±(r,t) and G»(r, t) are roughly similar to the shapes of 
these functions in the theory of homogeneous turbulence, as illustrated in Fig. 3.2 of chapter 
3 in Ref. For a more detailed description of these functions, we refer to J3T]]. 



We now return to the structure factors, and compare Eq. (|5^) with results from a two- 
dimensional molecular dynamics simulation of iV = 50000 inelastic hard disks at a low area 
fraction = ^nricr 2 = 0.05, and a coefficient of normal restitution a = 0.85. Fig. 2a shows 
the prediction of Eq. (|52"|) in the low density limit at r = 30 collisions per particle, together 
with simulation results for S±(k,t) and S\\(k,t), as obtained by performing a Fourier trans- 
form on the momentum fields, coarse-grained into 256 x 256 boxes, followed by a circular 
average. At the corresponding time t, Eq. (p3|), gives the slightly smaller value r = 28.7 
for the number of collisions per particles. Already at this density, Enskog theory gives a 
quantitatively significant increase of the collision frequency by a factor \ ~ 1.08, leading to a 
predicted number of collisions per particle r = 29.7, for the case under consideration. For a 
more detailed comparison of both S±(k,t) and S»(k, t) at general densities, we refer to |]3T 
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Since the structure factor of transverse velocity fluctuations, S±(k,t), shows structure at 
wavenumbers k < ~ 0.06cr _1 in the ring kinetic theory, the spatial velocity correlation 
functions G±(r, t) and G»(r, t) (shown in Fig. 2b) also show structure up to and beyond dis- 
tances of the order 2tt£ ~ 106a, which is large compared to the Boltzmann mean free path 
/ ~ 6.3cr, but still small compared to the system size L = 886er. Furthermore, the simulation 
results show that the incompressibility assumption indeed holds in a range of wavenumbers 
V£|| < k < 1/^0) where the existence of a large distance cut-off length 27r£|| ^> 27r£ has been 
discussed in Ref. |31 |. The spatial velocity correlations measured in the simulations agree 
well with the calculated G±(r, t) and G\\(r,t) in the range 2tcIo < r < 2tt£\\, and exhibit an 
observable long range r~ d -tai\ for distances 27r^ < r < 27r^ii. At wavenumbers of the order 
of the minimal accessible wavenumber /c min = 2tt / L ~ 0.007cr _1 , effects from the periodic 
boundaries become important. To conclude, we observe that the more fundamental ring 
kinetic theory yields results identical to the more mesoscopic theory of fluctuating hydro- 
dynamics in their common region of validity, thus supporting the phenomenological theory 
presented in Ref. [E3ipT . 
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FIGURES 
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FIG. 1. Direct and restituting inelastic collisions with refection law vj 2 ■ <r = — avi2 • <x with 
< a < 1. 



(a) (b) 




FIG. 2. (a) S±(k,t) (filled circles) and S»(k,t) (open circles) as measured from molecular 
dynamics simulations of N = 50000 particles at an area fraction (j) = 0.05, coefficient of normal 
restitution a = 0.85 and number of collisions per particle r = 30, with Tq/ui = 1, compared 
with the prediction Eq. ( |52D from ring kinetic theory for S±(k,t) (solid line), (b) Corresponding 
measured Gn(r,t) and G±(r,t) (dashed lines), compared with the prediction (solid lines) from ring 
kinetic theory in the incompressible limit. 
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